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Section  I 

PLASMA  THEORY  AND  SIMULATION 


A.  Lower  Hybrid  Drift  Instability  Simulations  Using  ESI  Hybrid  Code 
Yu-Jiuan  Chen  (Prof.  C.  K.  Birdsall,  B.  I.  Cohen) 

The  saturated  electric  field  energy  spectrum  due  to  ion 
trapping  was  obtained  in  theory  ana  compared  with  simulation.  These  and 
all  of  the  simulations  will  be  given  in  an  ERL  report  now  in  preparation , 
to  be  submitted  to  the  Physics  of  Fluids. 

B.  Beaming  instabilities;  Magnetized  Rings 

Jin  Soo  Kim  (Prof.  C.  K.  Birdsall,  B.  I.  Cohen) 

Our  multi-ring  instability  study  (the  distribution  of  plasma  in 
velocity  space  perpendicular  to  the  magnetic  field  is  a  set  of  discrete 
rings)  was  begun  with  one,  two  and  four  rings.  The  one  ring  model  produces 
the  well  known  Dory-Guest-Harr i s  instability.1  The  multi-ring  model  corres¬ 
ponds  to  multi-velocity  neutral  injection,  as  well  as  to  initial  conditions 
common  to  simulation  models.  The  model  has  only  ion  rings  because  the  elec¬ 
tron  Larmor  radius  (in  terms  of  k^a^)  is  negligible  compared  with  that  of 
ions  (in  terms  of  k^.^1).  In  our  current  work,  the  rings  in  velocity 
space  have  equal  weight  (same  to*)  and  are  spaced  at  equal  intervals,  as 
shown  in  Fig.  t;  the  plasma  is  uniform  in  x  space,  as  is  the  magnetic  field. 

The  distribution  function  for  N  rings  is  the  sum  of  6-functions; 

fo(Vvl}  "  *£  2^6(vi‘v±s)  6(vI)  •  0) 

The  dispersion  relation  of  a  single  ring  for  perpendicular  propagation  (k|«0)  is 
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st ri but  Ion  of  ions  in  v  space.  The  ring  speeds  are: 

>r  2  rings,  Vx1  *  1-ot*  vi2  "  1+°5  ^or  **  rings,  ■  1-3a, 
2*1-0,  v^j*1+a,  v^*1+3a.  The  a  values  used  so  far 
e  0.00  to  0.06. 
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where  y.  •k.v./w  and  J  is  the  Bessel  function  of  the  first  kind  of  order  n; 
l  i  x  c  n 

2 

this  is  Eq.  (19)  of  Tataronis  and  Crawford  .  Hence,  for  the  multi-ring  model 
of  equal  density  rings  the  dispersion  relation  is  (ignoring  electrons) 


D(«fkx) 


i  H  *  1  d*r(u  )  nw 

p  y*  y  i  n  is  c 

M2"  "  y 

0)  s*1  n— «  is  dy,  a)  -  ntu 

c  is  c 


The  roots  of  DCco.k^)  «0  are  obtained  numerically,  using  the  code,  ROOTS  . 

Figures  2  and  3  shows  the  maximum  growth  rate  for  two  rings; 

2  2 

Figs.  A  and  5,  for  four  rings.  When  o»^/o>c  is  larger  than  the  threshold 

values,  some  of  the  Bernstein  modes  couple  to  each  other  and  the  waves  grow 

(y  =  Imag  (u)  ^0)  as  is  well  known,  y  Is  the  maximum  value  of  y  of  all 

max 

Bernstein  modes.  The  growth  rates  become  smaller  as  number  of  rings  goes 
from  1  to  2  and  to  4  by  spreading  the  distribution  over  v^ -space,  corres¬ 
ponding  to  warming  up  the  cold  ring,  a  result  found  earlier^.  As  is  shown 

2  2 

in  Fig.  3  and  Fig.  5  the  instability  threshold  values  of  u>  /u>  become  larger 

P  c 

as  the  number  of  rings  increases  from  1  to  2  to  A. 

These  studies  will  be  extended  to  larger  values  of  the  para¬ 
meters  and  then  to  non-uniform  envelopes  of  the  rings,  such  as  Maxwellian. 
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FIG.  5  Same  as  Fig.  3,  for  4  rings.  The  marginal  stability  values  of 

(a>  /to  )?  are  9*4  (a*0.04)  and  1 6 . 3  (a*0.06). 
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3.  M.  J.  Gerver,  "ROOTS,  A  Dispersion  Equation  Solver",  University  of  Cali¬ 
fornia,  Berkeley,  ERL  Memorandum  No.  M77/27  (1976). 

4.  M.  J.  Gerver,  C.  K.  Birdsatl,  A.  B.  Langdon  and  D.  Fuss,  "Normal  Modes 
of  a  Loss  Cone  Plasma  Slab  with  Steep  Density  Gradient",  Phy.  FI.  20, 
pp.  291-300,  Feb.  1977. 

C.  Field-Reversed  Plasma  Simulations,  Quas i neutral  ,  in  2d 
Douglas  Harned  (A.  Friedman,  C.  K.  Birdsall) 

A  realistic  simulation  of  field-reversed  layer  stability  problems 
requires  a  bounded  system,  normally  bounded  by  conducting  walls.  Addition¬ 
ally,  for  a  fusion  reactor,  the  region  near  the  walls  should  be  a  vacuum  or 
cold  gas.  The  replacement  of  the  vacuum  or  cold  gas  region  with  a  low  density 
plasma  is  inadequate  in  a  quasineutral  model.  When  such  a  replacement  is 
made,  large  electric  field  fluctuations  occur,  and  a  highly  restrictive  Courant 
condition  (At<Ax/v^,  where  v^  is  the  Alfven  speed)  is  imposed  in  the  low  den¬ 
sity  regions.  Our  one  dimensional  simulations  have  shown  that  large  amplitude 
waves  occuring  in  the  low  density  regions  can  distort  the  physics  in  the  dense 
part  of  the  plasma.  It  is  not  surprising  to  find  such  problems,  since  the 
quas i neutrali ty  assumption  becomes  invalid  for  low  density.  These  difficul¬ 
ties  seem  to  be  inherent  in  explicit  schemes  for  advancing  the  field  quantities. 
A  method  which  avoids  such  problems  in  low  density  regions  has  been  developed 
by  Hewett.^  Hewett's  method  implicitly  advances  the  fields  and  threats  the 
vacuum,  or  near  vacuum,  by  setting  the  resistivity  to  a  large  value. 

Boundary  conditions  impose  a  further  difficulty.  This  is  because 
a  quasineutral  code  does  not  model  the  details  of  the  sheath  region  at  the 
wall.  This  defect  prevents  the  specification  of  reflecting  wall  boundaries 
in  conjunction  with  conducting  wall  boundaries.  At  a  conducting  wall  we  know 


n  *  E 


=*  0  . 


(D 


While  this  condition  is  insufficient  to  solve  the  field  equat ions,  the  boundary 
condi t ion 

n  •  B  =  0  (2) 

implies  that  the  longitudinal  (curl-free)  and  transverse  (divergence-free) 
parts  of  the  tangential  electric  field  (E  and  E  )  must  each  vanish  at  the 

&  L 

conductor.  In  terms  of  the  potentials  A  and  <fr,  in  the  Coulomb  gauge,  we 
have  at  the  walls 

n  x  A  =  0  (3a) 

<}>  =  constant  (3b) 

% 

where  A  and  <p  are  such  that  B  =  V  x  A  and  E^  =  -V<p.  With  the  gauge  condition, 

V*A  =  0,  the  boundary  conditions  are  sufficient  to  determine  the  advance  of 
the  field  equations.  It  should  be  noted  that  in  one-dimensional  codes  (such 
as  our  QUADl)  the  boundary  problem  is  trivial  as  the  longitudinal  and  trans¬ 
verse  parts  of  the  electric  field  are  geometrically  decoupled. 

The  need  for  a  bounded  system  and  proper  treatment  of  low  density 
regions  have  motivated  us  to  change  the  field  solver  in  our  two-dimensional, 
doubly  periodic  code.  The  following  is  a  description  of  a  method  which  may 
handle  both  problems. 

The  quasineutral  field  equations  are 


E 


T-J—  (VxB)xB  -  —  JxB  +  t2-  VxB 
Airne  “  -  nec  •  •  hit 


(4) 
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where  n  is  the  resistivity,  J  the  ion  current  density,  and  the  electrons 
have  been  assumed  to  be  cold.  In  terms  of  the  vector  potential  A,  we  have 


E 


-1 

4irne 


Ax (vxA  + B  )  - 
-  -  -o 


—  Jx[(VxA)xB  ]  - 
nec  -  L  -  -o 


cn 

*7 


v2a 


(6) 


3A 

Tt 


(7) 


where  B  =  VxA+Bq.  B^  represents  a  constant  background  magnetic  field.  Equa¬ 
tions  (6  and  7)  can  be  combined  to  give 


3A 

n 


V2Ax(VxA+B  ) 

-  -  -o 


2  . 
4tt 


•  (8) 

t 


The  right  hand  side  of  Eq.  (8)  requires  a  decomposition  into  longitudinal 
and  transverse  parts.  In  special  cases  (e.g.,  one-dimensional  cases  or 
Hewett's  axisymmetric  model)  this  can  be  done  geometrically.  However,  in 
general,  and  for  the  non-ax i symmetri c  long  layer  problems  that  we  wish  to 
examine,  geometric  decomposition  is  not  possible. 

Equation  (8)  can  be  written  as 

“  ' c7*  ■  W7  [<VXA)*B01  (9) 

d  t 


where  <J>  represents  a  potential  such  that  This  equation  now  has  a 

form  somewhat  similar  to  that  of  the  Navier-Stokes  equation  for  incompres- 


A  critical  difference  between  the  two  equations  lies  in  the  boundary  condi¬ 
tions.  In  Eq.  (9),  the  potential  <J>  and  one  component  of  the  vector  potential, 
A,  are  known  at  the  boundaries.  In  the  Navi er-Stokes  equation,  the  pressure 
at  the  boundary  is  generally  unknown,  while  both  components  of  the  flow  velo¬ 
city,  v,  are  specified.  Therefore,  some  methods  used  in  solving  the  two- 
dimensional  Navier-Stokes  equation  ’  may  be  applied  to  Eq.  (10),  but  not 
without  significant  modifications  to  handle  the  different  boundary  condi¬ 
tions.  The  fact  that  the  scalar  potential  is  specified  on  the  boundary  in 
the  quasineutral  equation  makes  it  a  straightforward  process  to  obtain  <f> 
from  a  Poisson  equation  (as  opposed  to  the  fluid  case  where  the  boundary  con¬ 
dition  must  be  predicted  to  solve  a  Poisson  equation  for  P) .  Writing  the 
right-hand  side  of  Eq.  (9)  as  ^(A)  ,  we  have,  in  two  dimensions, 


3A 

if  -  cv  - 

(11a) 

3A 

-r*-  -  C7  <f>  *  <3T  (A) 

at  yr  ^y  - 

(11b) 

V2*  -  -v.g(A)  . 

(11c) 

Equations  (I1a,b,c)  form  a  non-linear  system  of  three  equations  in  three  un¬ 
knowns,  which  must  be  solved  simultaneously.  Equations  (11a  and  11b)  are  non¬ 
linear  in  that  they  have  products  of  A  and  A  ,  but  neither  contains  a  term 

-x  -y 

2  2 

like  A  or  A  .  These  equations  may  each  be  advanced  in  a  Crank-Nicolson 
x  y 

scheme.  This  forms  a  predictor-corrector  method  for  the  advance  of  all  three 
quantities.  For  the  m-th  of  M  iterations  we  have 
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An+1'm-An 

X  X 


-  "  f  (Gx«frn+1*m  +  G*")  +  j  <g"+1’m(A)+gr(A))  <12a) 


An+1 ,m_An 


|  (  Gy<(>n+1  ’m  +  Gv<frn  )  +  J  ^+1’m(A)  "^(A))  (12b) 


2  'vJy 


>n+1’m  -  -  |0h+1*"»(A)  . 


(12c) 


L,  D,  and  6  are  the  finite  difference  approximations  for  the  Laplacian,  diver¬ 
gence,  and  curl,  respectively.  The  operator  L  should  be  such  that 


L$  *  D(i<t>  . 


With  interlaced  grids,  using  four-point  operators  to  represent  the  deriva¬ 
tives  in  the  gradient  operators,  i.e.. 


i+U+1  vi+1,j-1  i”1 ,  i+1 


i+i» j+i 


)  .  (14) 


Eq.  (13)  implies  that  L  be  defined  by 


L*  5 


•+1,j+1 


i-H.j-1  *  h-1,J+1  * 
2 (Ax) 2 


Because  all  three  of  Eqs.  (12)  involve  the  solution  of  a  large  matrix  equation, 
the  method  is  totally  impractical  for  large  M.  It  is  assumed  that  M  < 2  wl 1  be 
sufficient.  While  the  time  requirements  for  this  type  of  solution  might  be 
prohibitive  in  a  fluid  code,  the  largest  fraction  of  time  used  by  a  hybrid 
code  is  expected  to  be  in  the  particle  mover.  Hopefully,  this  will  allow  the 
flexibility  to  solve  a  system  of  equations  such  as  (12).  If  a  technique  such 
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as  orbit  averaging  is  used  for  the  particle  mover,  then  the  time  constraints 
on  the  field  solver  will  be  reduced  further. 


The  most  attractive  way  to  solve  Eqs.  (12)  appears  to  be  the 

ICCG  (incomplete  Cholesky  conjugate  gradient)  method,  because  of  the  large 

variation  in  the  magnitude  of  the  diagonal  elements  from  the  plasma  to  the 

4 

vacuum  regions.  The  first  predictor  steps  should  require  a  relatively 
small  number  of  iterations,  as  it  would  be  unnecessary  to  demand  a  small 
residual  in  the  predictor  solution. 

Once  A  and  <f>  have  been  advanced,  the  new  values  for  E  and  B 
can  be  obtained  from 


B 

E 


n+1 

n+1 


VxAn+1  +  B 
-o 

igr<A)"+' . 


These  fields  can  then  be  used  to  advance  the  particles  with  the  standard 
techniques. 
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Effects  of  Intrinsic  Orbital  Stochasticlty  on  Resonant 
Microinstability  (Bruce  Cohen  and  Alex  Friedman) 


Complex  inhomogeneous  equilibria  can  exhibit  Intrinsic  stochasticlty; 
i.e.  certain  classes  of  particles  have  stochastic  orbits.  An  illustrative 
and  familiar  example  arises  In  simple  mirror  systems.  For  axlsymmetrlc 
mirror  machines  particle  energy  &  and  canonical  angular  momentum  Pg  are 
conserved  in  the  absence  of  turbulence  or  other  perturbations  that  break 
azimuthal  symmetry  and  Introduce  explicit  time  dependence.  When  the  Larmor 
radius  is  small  compared  to  the  magnetic  scale  lengths,  the  magnetic 
moment  u  and  the  longitudinal  action  J  are  approximate  (or  adiabatic) 
Invariants.  The  magnetic  moment  can  experience  changes  when  somewhere  along 
the  particle  trajectory  the  cyclotron  frequency  or  Its  harmonic  equals  a 
harmonic  of  the  axial  bounce  frequency.  For  changes  In  u  from  one  bounce  to 
another  that  are  unrelated,  the  particle  motion  remains  orderly  and 
deterministic,  but  the  changes  In  u  as  the  particle  passes  through  resonance 
on  successive  axial  transits  can  be  viewed  as  a  stochastic  process. 
Neighboring  stochastic  orbits  in  the  phase  space  of  y,  J,  and  their 
conjugate  angles  appear  to  diverge  exponentially  in  time  when  averaged  over 
many  bounce  periods. 

What  effects  does  intrinsic  orbital  stochasticlty  have  on  micro¬ 
instability?  Intrinsic  stochasticlty  is  not  collective  In  nature;  the  orbit 
separation  of  various  sets  of  neighboring  particles  in  phase  space  is 
uncorrelated  when  viewed  over  many  bounce  periods.  Furthermore, 
stochasticity  does  not  influence  the  moments  of  the  distribution  function 
and,  in  particular,  does  not  alter  any  free  energy  possibly  available  nor 
the  macroscopic  charge  densities  and  currents.  Therefore,  we  conjecture 
that  the  effects  of  stochasticlty  on  mlcrolnstablllty  may  be  weak  provided 
that  the  equilibrium  orbits  are  not  too  much  distorted.  The  following 
discussion  fills  in  the  details  of  this  argument  and  supports  Its  conclusion. 

For  simplicity  we  assume  that  the  linear  growth  rate  of  an  unstable 
mode  and  the  stochasticlty  rate  (orbital  separation  rate)  are  both  much 
smaller  than  the  mode  frequency,  cyclotron  and  bounce  frequencies.  Upon 
integrating  the  Vlasov  equation  along  its  characteristics,  I.e.  performing 
the  Integral 

3f  •*  /  dt '  exp  Clky(t')  -  lut'] 

•/  -00 
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(1) 


in  slab  coordinates,  denominators  appear  with  resonance  condition 


Aw  =  w  -  £w  ci(J,u,X)  +  2p<^(u)  +  kVd(u)  -  0  ,  (2) 

where  is  the  average  cyclotron  frequency,  is  the  bounce  frequency, 

Vd  is  the  drift  velocity,  and  X  is  the  guiding  center  position  (see  Smith, 
Byers,  and  LoDestro  UCRL-82674). 

We  now  presume  that  the  adiabatic  invariants  u  and  J  begin  to  slowly 
diverge  due  to  intrinsic  stochasticity.  Hence,  Aw  *  Aw(t)  and  a  particle 
will  cease  to  be  resonant  when 

f  S dt '  Aw(  t ' )  ■  9(  tt)  , 

JQ 

which  defines  and  relates  it  to  the  stochasticity  rate.  Over  the  time 
period  the  particle  can  do  work  on  the  wave  and  vice  versa.  Meanwhile, 
however,  accompanying  the  loss  of  resonance  for  one  particle,  there  Is  an 
equal  probability  that  a  neighboring  particle  in  phase  space  is  coming  into 
resonance.  Resonance  again  persists  over  a  time  period  of  order  t$.  The 
1  inear  perturbations  of  the  resonant  particles  will  add,  and  the  linear 
dielectric  response  will  be  largely  unaffected  by  stochasticity.  The  crux 
of  the  argument  is  that  with  stochastic  orbits  we  expect  an  equal  flux  of 
particles  into  and  out  of  resonance  at  any  point  on  the  separatrix  between 
resonant  and  nonresonant  particles  but  the  wave-particle  interaction  while 
in  resonance  is  unaffected.  Nonlinear  aspects  of  the  wave-particle  inter¬ 
action  are  decidedly  affected  when  intrinsic  stochasticity  limits  the 
duration  of  resonance,  essentially  because  nonlinear  effects  associated  with 
particles  coming  in  and  out  of  resonance  are  not  simply  additive. 


Cartoon  of  separatrix  for  particles 
resonant  with  an  Ion  cyclotron  flute 
mode  (see  Smith  et.  al.,  UCRL-82674) 
showing  fluxes  of  stochastic 
particles. 
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When  x$  becomes  so  short  that  x"1  is  comparable  to  the  smallest 
characteristic  frequency  helping  to  determine  the  resonance  condition, 

Eq.  (2),  then  the  resonance  is  destroyed.  In  this  case  the  particle  will 
interact  simultaneously  with  many  of  the  bounce  harmonics,  and  there  is  no 
resonance  with  any  one  particular  bounce  harmonic.  Thus,  stochasticity  will 
profoundly  influence  the  linear  aspects  of  resonant  microinstability  when 
the  stochasticity  rates  becomes  comparable  to  the  characteristic  frequencies. 

Random  collisions  have  a  very  different  effect  on  microlnstability. 

Use  of  the  simplest  Krook  model  yields  the  following  modification  to  the 
resonant  denominators  appearing  in  6f, 

b  siii  .  +  2pu^  +  kVd  +  iv  .  (3) 

* 

In  the  absence  of  collisions,  Im  Aw  3  Imw  3  y,  the  linear  growth  rate.  With 
collisions,  Im  Aw  »  y  +  v  .  Thus,  collisions  will  have  a  significant 
effect  when  »  9^y).  Physically,  collisions  produce  a  real  diffusion  of 
particles  in  phase  space  and  destroys  and  P0  conservation,  as  well  as  alter 
u  and  J.  Wave  growth  is  necessarily  accompanied  by  momentum  and  energy 
exchange  with  resonant  particles.  Therefore,  at  the  point  when  the  collision 
rate  becomes  large  enough  to  influence  the  energy  and  momentum  exchange  of 
the  resonant  particles  interacting  with  the  wave,  microinstability  would  be 
significantly  altered. 

To  illustrate  collisional  effects,  consider  the  force  on  a  resonant 
particle  due  to  a  wave  in  the  presence  of  randomizing  collisions. 


at  «*  *  s  e 


-lAwt 


+  c.c.  -v  6v 
c 


where  6v  is  the  perturbed  velocity  in  the  frame  of  the  unperturbed  particle 
motion,  F  Is  the  amplitude  of  the  force,  and  Aw  might  be  given  by  Aw  3  w-kv, 
Eq.  (2),  or  one's  favorite  resonance  condition.  The  solution  of  Eq.  (4)  is 


_  -1  -iAwt 

6v  3  F  e  +  c.c. 
-1Aw+\£ 


-  i  Awt  +  1  <$ 


a +(VC+Y  Y 


+  c .  c .  , 


where  cos  6  *  (vc  +y)/^  +  (vc  +  y)^»  A  s  Re  Aoj,  and  any  initial 
transient  has  been  allowed  to  damp  away.  It  is  obvious  from  this  simple 
expression  that  the  energetics  of  the  wave-particle  Interaction  are 
significantly  affected  when  v.becomes  comparable  to  y. 

It  appears  that  there  can  be  considerable  stochasticlty  without  there 
being  much  effect  on  the  linear  aspects  of  resonant  microinstability.  In 
simulations  such  as  Friedman's  three-dimensional  linearized  simulations  of 
field-reversing  rings  and  mirrors,  intrinsic  stochasticlty  of  particle 
orbits  is  often  observed.  Based  on  the  arguments  presented  here.  If 
particle  statistics  in  the  simulation  are  very  good,  but  not  economically 
unmanageable,  so  that  there  is  the  necessary  cancellation  of  the  stochastic 
particle  flux  into  and  out  of  resonance  with  the  wave,  then  one  expects  that 
Friedman's  code  could  investigate  resonantly  driven  instabilities  even  In 
the  low  growth  rate  regime.  Obviously,  the  higher  the  growth  rates  of  the 
modes  under  investigation  the  less  perfect  the  cancellation  of  stochastic 
particle  flux  need  be  and  the  fewer  simulation  particles  necessary. 

Remaining  to  be  done  is  the  quantitative  assessment  of  the  particle 
statistics  requirements  for  simulation  of  modes  in  the  regime  where  the 
stochasticity  rate  exceeds  the  instability  growth  rate. 

We  thank  John  Finn  and  Jim  Albritton  for  illuminating  discussions  of 
these  issues,  and  Brendan  McNamara  for  posing  the  problem. 
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E.  Effects  of  Numerical  Dissipation  on  Simulation  of  Fast  and  Slow 
Space-Charge  Waves 

Bruce  I .  Cohen 


The  introduction  of  backward  time-differencing  in  the  numerical 
solution  of  either  ordinary  or  partial  differential  equations  with  time 
dependence  can  lead  to  frequency-dependent  numerical  dissipation.  This 
can  be  exploited  and  provide  a  means  of  performing  digital  filtering  in 
numerical  simulations  of  time-dependent  phenomena  [1-3].  In  Ref.  1 
digital  time  filtering  is  combined  with  an  Implicit  differencing  scheme  to 
achieve  damping  of  high  frequency  waves  in  an  unconditionally  stable 
algorithm  that  allows  use  of  large  time  steps.  This  approach  is  suitable 
for  solution  of  fluid  equations  and  linearized  kinetic  equations,  e.g., 
the  drift-kinetic  equation  [1].  Unfortunately,  implicit  solution  of 
nonlinear  kinetic  equations  frequently  leads  to  inversion  of  large 
non-sparse  matrices^]. This  is  particularly  true  of  particle  codes,  and 
thus  the  implicit  methods  described  in  Ref.  1  have  not  been  applied  to 
particle  simulations.  Nevertheless,  digital  filtering  has  been  frequently 
used  in  particle  codes  to  damp  unwanted  high  frequency  modes  and  recently 
has  been  combined  with  a  new  technique  called  orbit-averaging  to  allow  the 
use  of  a  large  time-step  in  the  solution  of  the  self-consistent  fields  [2]. 
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An  interesting  and  important  question  has  been  raised  regarding  the 


use  of  numerical  algorithms  with  artificial  dissipation  added.  Because 
negative  energy  modes  in  plasmas  or  in  other  dielectric  media  can  be 
destabilized  by  dissipation,  it  is  wondered  whether  numerical  dissipation 
can  artificially  destabilize  negative  energy  modes.  The  following 
discussion  demonstrates  that  numerical  dissipation  can  in  fact  destabilize 
negative  energy  modes,  but  that  with  some  ingenuity  dissipative  algorithms 
can  be  devised  that  damp  both  positive  and  negative  energy  modes 
increasingly  as  the  mode  frequency  increases. 

The  simplest  model  admitting  positive  and  negative  energy  modes  is 
that  of  a  cold,  unmagnetized,  drifting  plasma.  The  infinite  homogeneous  ' 
plasma  dispersion  relation  for  plasma  waves  becomes 

(w  -  kuQ)2  *  up  ,  (la) 

with  solutions 

m  3  kug  ±  Up  ,  (lb) 

where  Uq  is  electron  drift  speed  relative  to  an  immobile  neutralizing 
ion  background,  k  is  the  wavenumber,  and  Wp  is  the  electron  plasma 
frequency.  The  +(-)  sign  in  Eq.  (lb)  corresponds  to  the  fast  (slow)  beam 
mode,  which  has  positive  (negative)  wave  energy.  As  we  shall  illustrate, 
the  slow  wave  can  be  driven  unstable  by  dissipation. 

Consider  the  one-dimensional  particle  simulation  algorithm  given  by 


vn+l/2  ,  vn-l/2  _  e£nAt/m 

(2a) 

xn+l  ,  xn  ■»*  Vn+1/,2At 

(2b) 

(1-e)  En  +  e|j«  En-1  ■  -  4ire(nn  -  nQ) 

(2c) 

>  - 
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where  X  and  V  are  electron  positions  and  velocities,  E  is  the  electric 
field,  ng  is  the  unperturbed  homogeneous  background  number  density,  the 
superscript  n  denotes  the  time  level,  and  e  is  a  centering  parameter  (e  *  0 
is  perfectly  time-centered).  To  analyze  the  system  of  equations  given  in 
Eq.  (2),  we  linearize,  use  the  relation 


n^  *  -  ng3X^/3x 


employ  the  Fourier  representation 


E  *  E  exp(-  i cot  +  ikx)  +  c.c. 


for  all  perturbed  quantities,  and  define 


A  =  exp(-  iuAt/2  +  ikUgAt/2) 


to  obtain 


(A  -  A_i)V  =  -  eAtE/m 


(A  -  A'1)?  =  VAt 


(1  -  e  +  eeia)/it)E  =  4TrennX 


The  dispersion  relation  for  this  system  of  equations  is  easily  determined: 


(1  -  e  +  eeia)At)(A  -  A"1)2  +  u>2At2  =  0 


For  |wAt|  «  1,  e1loAt  «  1  +  iwAt  and  Eq.  (5)  has  solutic 


A2  a  1  ±  iot  -  §j—  ±  +  ...  , 
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where 


2,2^2, 


a  =  aJpAt/(  1  +  e  w  At  )  «  1 


From  Eq.  (6)  and  the  definition  of  X,  we  deduce 


Re  uj  «  kuQ  +  oip/(  1  +  e2w2At2) 


I*  I*  1 1 


eaxjjpAt 

2(1  +  e2(*£  At2) 


One  concludes  from  Eqs.  (7a)  and  (7b)  that  the  fast  mode  is  damped 
and  the  slow  wave  destabilized  for  £  >  0  and  kug  ±  ^  >  0.  Furthermore, 

| Im  cd |  «  ( l/2)eoipAt | Re  uj|,  i.e.,  the  damping  and  growth  rates  are  frequency 
dependent.  The  same  sorts  of  results  are  obtained  with  use  of 


ir  *  -  <’  ♦  «>  4"e(%  -  "o>  -e  h  e"'1 


in  place  of  Eq.  (2c). 

A  useful  interpretation  of  the  effects  of  finite  e  in  Eqs.  (2c)  and  (8) 
can  be  given  by  observing  that  these  equations  can  be  cast  in  a  form 
equivalent  to 


§£  (ff  +  4ttJ  +  47raE)  =  0 


where  J  is  the  plasma  current  and  a  is  an  electrical  conductivity 
representing  the  interaction  of  the  plasma  with  a  resistive  background. 
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For  | a/a) |  «  1,  the  dispersion  relation  given  in  Eq.  (1)  becomes 


=  kug  +  - 


1  + 


(  4tto  )  ‘ 


(ku0  +  wp)‘ 


1/2 


iTT2g  \ 

ku0  -  WP  / 


(10) 


For  a  >  0  and  kug  >  up,  the  slow  wave  is  unstable  and  the  fast  wave  is 
damped.  This  is  the  underlying  physical  mechanism  for  the  resistive  wall 
ampl ifiert?].  The  structure  of  Eq.  (10)  is  very  similar  to  Eqs.  (7a)  and 
(7b)  and  permits  an  identification  of  a  frequency-dependent  numerical 
conductivity. 

An  alternative  algorithm  to  that  given  in  Eq.  (2)  is 

(1  -  e)(Vn+1/2  -  Vn'1/2)  e(Vn+1/2  -  Vn"3/2)  =  -eEnAt/m  (11a) 

Xn  -  Xn-1  =  Vn_1/2At  (lib) 


Jj=-  =  -  4ire(nn  -  nQ)  .  (11c) 

With  e  >  0,  the  time  derivative  of  the  velocity  is  backwards-differenced 
with  respect  to  the  electric  field  in  Eq.  (11a).  Employing  the  definitions 
and  Fourier  representation  introduced  earlier,  we  linearize  and  obtain  the 
bi-cubic  dispersion  relation 

\Z{\ 2  -  l)2  -  |  (X2  -  l)3  +  u)p2At2X4  =  0  .  (12) 

2 

For  WpAt,  «  1  Eq.  (12)  has  one  branch  with  solution  given  by  A  =  -  e/(2  -  e), 
which  is  damped  for  e  <  1.  For  u>pAt,  |e|  «  1,  the  remaining  solutions  are 

2  ..2 

o  ^  p  o 

r  »  1  +  oipAt  (1  +  ieojpAt/4)  -  -  (1  +  iaopAt/4)4  +  ...  , 
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or  equivalently 


2 

ui  a  kug  +  Wp  -  ieujp  At/4 


(13a) 


\X2\  «1  -  ewp2At2/4  <  1 


(13b) 


for  e  >  0.  Thus,  both  fast  and  slow  beam  modes  are  damped.  The  damping 

•p 

rate  is  eiopAt/4  and  vanishes  in  the  limit  that  e  ■*  0.  This  algorithm 
has  an  effective  numerical  conductivity  which  is  positive  for  the  fast  mode 
and  negative  for  the  slow  mode. 

For  copAt  >  >  1,  Eq.  (12)  becomes 

(1  -  e/2)x4  a  -  tup  W 

for  one  branch  of  solutions;  and  thus  (\|  a<DpAt/(l  -  e/2}^2  >  1 
for  e  <  2,  i.e.,  the  algorithm  is  numerically  unstable  for  large  time 
steps.  The  remaining  solutions  of  Eq.  (12)  are  given  by 


,2  -  _J _ .  1/  1  .  2s 

2<op2At2  '  2  \*.p4At4  "P2«2 


for  <^pAt2  »  e,  1  and  are  heavily  damped  (|A2|«  1). 

For  e  *  0  and  general  value  of  ^pAt,  Eq.  (12)  simplifies  to  the 
dispersion  relation  for  the  conventional  leap-frog  algorithm. 


(A2  -  l)2  +  wp2At2A2  =  0 


with  solution 


sin2(kug  -  w)t  =  ^p2At2/4  , 


which  is  stable  for  WpAt  <_  2.  For  small  but  finite  e  and  a  =  oipAt  very 
nearly  equal  to  2,  the  dispersion  relation  in  Eq.  (12)  can  be  expanded  for 
small  perturbations  <5e  and  6a,  with  respect  to  e  =  0  and  a  =  2, 

D(A  ;  e,  a)  555  D(A.  =  -  1;  e  =  0,  a  =  2)  +|5e'|“j  +  =  0  ,  (14) 

where  D(*2;  e,  a)  is  given  by  Eq.  (12)  and  D  =  0  at  X2  =  -  1  for  e  -  0 
and  a  =  2.  Equation  (14)  then  yields 


This  demonstrates  that  the  introduction  of  a  small  amount  of  biasing, 

0  <  6e  «  l,  shifts  the  stability  boundary  of  the  algorithm  to  slightly 
smal ler  values  of  wpAt. 

A.  B.  Langdon  has  contributed  the  valuable  observation  that  the 
backwards  biasing  in  Eq.  (11a)  occurs  in  a  Lagrangian  equation  and  is 
invariant  under  a  Galilean  transformation.  Thus,  dissipation  is  introduced 
in  a  way  that  is  independent  of  reference  frame;  and  artificial 
destabilization  of  a  beam  mode  should  not  be  expected  when  a  drift  is  added. 
In  contrast,  the  backwards  biasing  in  the  first  algorithm,  Eq.  (2),  is 
performed  in  the  field  equation,  which  is  Eulerian  in  nature.  The  resulting 
effects  are  not  invariant  under  a1  Galilean  transformation,  and  the  mischief 
caused  in  the  slow  beam  mode  is  not  so  surprising. 

Birdsall  has  commented  that  biasing  the  field  equation  in  Eq.  (2c)  is 

2 

equivalent  to  changing  the  dielectric  medium  in  such  a  way  that  Up  becomes 
complex.  This  results  in  the  destabilization  described  in  Eq.  (10). ^ 
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However,  with  the  biasing  indicated  in  Eq.  (11a),  the  equation  of  motion  is 
changed  in  a  purely  dissipative  way,  effectively 


dv  dv 
dt  "*  dt 


e' 


=  -  eE/m 
dt^ 


where  e'  =  eAt/2  >  0;  and  upon  Fourier  analyzing,  (w  -  kug)  -*■  £oj  -  kug  + 
ie"(u  -  kug)  ].  This  scheme  should  damp  all  waves  for  small  values  of  oipAt. 

One  concludes  that  digital  time  filtering  need  not  lead  to  artificial 
destabilization  of  negative  energy  waves.  Algorithms  can  be  devised  which 
damp  both  positive  and  negative  energy  waves  as  an  increasing  function  of 
mode  frequency  and  time  step.  However,  in  an  explicit  differencing  scheme, 
the  backwards  biasing  of  the  difference  equations  necessary  to  achieve  time 
filtering  is  likely  to  make  the  time-step  constraint  for  numerical  stability 
slightly  more  stringent. 

I  am  pleased  to  acknowledge  many  helpful  discussions  with  8ill  Fawley, 
Bruce  Langdon,  and  C.  K.  Birdsall,  and  am  indebted  to  Brendan  Godfrey  for 
posing  the  problem.  This  work  was  performed  under  the  auspices  of  the  U.S. 
Department  of  Energy  at  Lawrence  Livermore  Laboratory  under  contract  number 
W- 7405 -ENG-48. 
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F.  Orbit  Averaging:  A  Technique  for  Filtering  out  High  Frequency  Waves 
and  Fluctuations. 

Vincent  Thomas  (Prof.  C.  K.  Birdsall  and  B.  I.  Cohen) 

A  technique  first  used  by  Cohen  et  a?.  (Ref.  l),  referred  to  as 
orbit  averaging,  has  been  written  for  use  in  ESI,  our  1-D  electrostatic  code 
in  an  attempt  to  filter  out  high  frequency  electron  waves  and  leave  the  low 
frequency  motion  of  both  particle  electrons  and  ions  and  the  associated  self- 
consistent  electric  field.  In  this  method,  the  electrons  and  ions  are  ad¬ 
vanced  on  their  respective  time  scales,  but  the  fields  are  solved  for  only 
on  the  ion  time  scale.  The  field  solver  uses  a  time  filtered  source  term, 
the  charge  density,  in  order  to  filter  out  the  high  frequency  phenomena. 

For  all  cases  considered  here,  the  time  filtering  consists  of  an  equally 
weighted  average.  The  flow  chart  is  shown  in  Fig.  1.  Also  shown  is  the 

scheme  for  moving  from  one  macro  time  step  to  the  next  macro  time  step.  The 

velocities  of  the  electrons  are  known  at  the  short  hash  marks,  and  the  velo¬ 
cities  of  the  ions  are  known  at  the  long  hash  marks.  The  fields  are  solved 
for  at  the  dots.  After  the  fields  are  solved  for,  the  electrons  and  ions 

are  mowed  from  NAT-iAt.  to  (N+l)AT-iAt.  and  from  NAT-iAt  to  (N+l)AT-iAt  , 

ii  e  e 

respectively.  The  electrons  must  be  moved  in  a  series  of  short  steps  to 

satisfy  u)  At  < 2.  The  electron  time  step  and  the  ion  time  step  will  be  re- 

pe  e  r 

ferred  to  as  the  micro  step  and  the  macro  time  step,  respectively. 

In  attempting  to  introduce  some  flexibility  one  may  introduce 
various  biasing  parameters.  The  first  of  these  is  in  the  field  solver  where 
one  uses  Eq.  (l) 

-eV^<J>N+^  -  (1-e)V^<j>N  ■  47re(n.(x)  -  n^x))  (l) 
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initialization 


When  was  larger  than  u^,  simulations  of  a  cold  plasma 

produced  clear  cut  oscillations,  with  frequencies  at  the  lower  hybrid 

frequency.  Runs  were  made  for  oo  *5  and  2.5  with  u  *1,2,3,  and  4. 

ce  pe 

For  u>  less  than  u>  ,  the  oscillations  were  not  as  clear,  and  energy 
ce  pe 

conservation  was  poor. 

Figure  2  shows  an  m  vs  k  diagram  for  a  particular  set  of  para¬ 
meters.  The  formula  used  to  calculate  the  theory  curve  means  using  the 
usual  finite  kAx  correction  appropriate  to  linear  weighting, 


a)  (kAx) 


(2) 


The  difference  between  theory  and  simulation  is  only  significant  for  modes 
14  and  15.  No  k-space  smoothing  was  used  in  any  of  these  simulations  and 
so  some  of  the  simulations  were  of  poor  quality.  The  frequencies  were 
clearly  visible  in  the  electrostatic  energy  history  plots  in  all  cases,  as 
the  simulations  were  done  by  exciting  only  one  mode  at  a  time.  The  modes 
were  excited  with  a  sinusoidal  perturbation  in  ion  and  electron  positions. 

For  cold  lower  hybrid  waves,  it  was  necessary  to  use  e  larger 
than  1  to  achieve  satisfactory  energy  conservation.  Using  e  greater  than 
1  amounts  to  an  interpolation  of  the  potential  at  a  time  later  than  the 
(N+1 )  ^  macro  time  step.  For  0.6  <  e  <2.5  the  period  of  the  hybrid  osci  1  1a- 
tions  observed  was  within  a  few  percent  of  the  expected  period,  assuming 
the  usual  kAx  correction  due  to  linear  weighting.  Figure  3  shows  the  vari¬ 
ation  of  total  energy  vs.  e,  with  large  loss  for  small  e.  The  field  energy 
and  the  kinetic  energies  damp  at  approximately  the  same  rate;  this  means 
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to  solve  for  <j>N+\  Here  N+1  designates  the  (N+1)^1  macro  time  step  which, 
in  general,  could  be  different  from  the  ion  time  step;  for  all  work  done 
here  the  macro  time  step  is  equal  to  the  ion  time  step,  n. (x)  is  the  ion 
charge  density  at  macrotime  steps;  n^(x)  ■  $  the  electron  charge  density 
averaged  over  all  of  the  micro  time  steps  from  NAT-JrAte  to  (N+1)AT-iAte 
(see  Fig.  1 ) . 

Another  biasing  parameter  can  be  introduced  in  the  electron 
mover.  For  the  predictor  step,  the  electrons  are  moved  by  E®-?$N,  which 
is  taken  to  be  constant  in  time.  For  the  corrector  iteration,  one  could 
use  an  E(t)  which  is  linearly  interpolated  from  the  (N)^1  to  the  (N+1)^1 
macro  time  step,  as 

E (t)  *  -  (l-bias2)v<|>N  +  (t-NAT)  (1+bf as2)  ( 

where  AT  is  the  macro  time  step.  In  the  simulations  presented  in  this 
report  bias2=0.  The  effects  of  the  parameter  bias2  will  be  commented  on 
at  a  later  time. 

The  primary  advantage  of  orbit  averaging  is  the  reduction  in  the 
number  of  particles  needed  for  a  simulation  [1],  This  is  possible  because 
the  high  frequency  part  of  the  thermal  spectrum  is  removed  by  this  simula¬ 
tion  scheme. 

Simulations  have  been  made  for  cold  and  warm  magnetized  plasmas 
and  for  warm  unmagnetized  plasmas  using  real  mass  ratios  (e.g.,  mj/me*l836). 
Preliminary  results  have  been  encouraging.  Results  of  the  three  different 
cases  will  be  discussed  separately,  with  the  results  of  the  warm  unmagnetized 
simulations  being  deferred  until  a  later  date. 
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Frequencies  of  oscillations  for  a  cold  magnetized  plasma.  The 

parameters  were:  u>  *2.0,  u>  “5.0,  m./m  *  1760,  the  electron 
pe  ce  i  e 

time  step  was  Ate*.03,  the  ion  time  step  At.  *3*0,  there  are 
512  electrons  and  512  ions,  NG  ■  32,  and  the  amplitude  of  the 
sinusoidal  perturbation  is  'vIO  ^  of  a  particle  separation. 

The  bias  parameter  e  is  equal  to  1.  The  curve  represents  theory, 
the  points  simulation  data.  Linear  weighting  was  used. 


TE  (T) 

TeToI 


s 


i 


FIG.  3 


Dependence  of  energy  conservation  on  bias  parameter  e. 

For  these  simulations  At  -3.0,  At.  -3.0,  m./m  -  1836,  u>  -2.0, 

6  I  16  p6 

uce“5.0,  NG  -  32 ,  L-32,  and  there  are  512  cold  electrons  and 

ions.  The  amplitude  of  the  initial  sinusoidal  perturbation 

was  ^ 10  ^  of  a  particle  separation.  T.E.(T)  and  T.E.(O)  are 

the  total  energies  at  T(-<u  t-300)  and  at  t-0,  respectively. 

pe 
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that  the  low  frequency  energy  (e.g.,  field  energy  at  u)^u>IU)  damps  at  the 
same  rate  as  the  higher  frequency  energy  (e.g.,  electron  kinetic  energy). 
Representat i ve  plots  of  one  run  are  shown  in  Figs.  (4  and  5).  The  exchange 
of  electrostatic  energy  and  ion  kinetic  energy  is  very  clear,  as  is  repre¬ 
sentative  of  ion  plasma  oscillations.  When  the  plasma  is  made  warm,  lower 
hybrid  oscillations  are  again  observed.  In  addition,  far  better  energy 
conservation  was  obtained  (less  than  ]%  error  up  to  w^t  *  300)  for  all  e 
tried  from  .6  to  1.  At  present  we  have  not  included  spectral  analysis  on 
this  code  and  therefore  large  (and  strongly  nonlinear)  initial  excitation 
was  required  to  produce  measurable  oscillations. 

The  next  step  is  twofold:  first,  to  do  extensive  parameter 
studies  and  to  understand  the  effects  produced  by  those  parameters,  both 
physical  and  nonphysical;  secondly,  to  add  suitable  diagnostics  so  that 
linear  waves  in  warm  plasmas  may  be  studied.  More  physical  simulation 
models  could  also  be  used,  such  as  k*B^0. 

REFERENCES 

1.  B.  I.  Cohen,  T.  A.  Brengle,  D.  B.  Conley,  and  R.  P.  Freis,  HAn  Orbit- 
Averaged  Particle  Code11,  Lawrence  Livermore  Laboratory  Report,  UCRL- 
82832  (June  1979).  Accepted  for  publication  in  J.  Comp.  Phys. 
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Same  parameters  as  in  Fig.  4.  (a)  electron  kinetic  energy  vs 

time  from  t  ■  0  to  t *  300  (b)  ion  kinetic  energy  from  t » 0  to 

t  *  300  (c)  total  energy  from  t  *  0  to  t*  300. 


Station  II 

CODE  DEVELOPMENT  and  MAINTENANCE 


A.  ESI  Code 

See  Orbit  Averaging  in  Sec.  I,  Part  E. 

B.  EMI  Code 

No  special  progress  to  report. 

C.  EZOHAR  Code 

No  special  progress  to  report. 

D.  RINGHYBRID  Code 

This  code  is  now  being  run  with  i n i t i al i zat ion  provided  by 
RIGIDROTOR  code  as  described  in  the  section  following,  in  the  Sherwood 
abstract  later  in  this  QPR,  and  in  an  ERL  report  in  preparation. 
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E.  RIGIDROTOR  Code  (Field  Reversed  Equilibrium  Solver) 


Alex  Friedman 

The  following  is  a  copy  of  the  LIBRIS  abstract  describing  the 
RIGIDROTOR  code.:  a  brief  description  of  the  applications  of  this 
code  appears  elsewhere  in  this  QPR. 


77  Alex  Friedman  RIGIDRGTGR 

DATE-  march,  1930  04/08/30  08:44:01 

rigid  rotor  equilibrium  solver  and  particle  code  initialization  program 
j.  c.  berKeley 

electronics  research  laboratory,  cory  hall 
university  of  California,  berKeley  ca  94613 
415-642-3477 

cdc  7600  fortran  1000  1  ines  in  use 


AB5TRACT- 

"r ig idrotor "  is  a  pacKage  which  calculates  field  reversed  equilibria 
using  the  t ime- independent  vlasov  equation,  and  outputs  the  fields  and 
a  set  of  particle  initial  conditions  so  that  these  equilibria 
may  be  used  in  particle  simulations.  the  pacKage  was  designed  for  use 
with  the  "rirghybrid"  linearized  3d  stability  code,  but  is  suitable 
for  other  axisymmetric  particle  or  hybrid  codes  with  minor  modification. 

in  the  present  impl ementat ion,  exponential  rigid 
rotor  vlasov  equilibria,  with  distribution  functions  of  the  form 
f  *  c  exp  C  -  (  h  -  omega  *  ptheta  )  /  t  ],  are  computed  on  an  r — z  grid 
U3ing  an  iteration  scheme  similar  to  that  of  sparks  and  finn. 


V  - 
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the  equation  solved  is  : 


C  psi  -  thdot  bnot  rsquared  /  2  ) 


1  a  -  r  dens  thdot  bnct  exp  (  -  ) 

(  psi max  ctemp  ) 


where  a  =  vector  potential  a theta,  psi  *  r  a,  psi max  3  ps i (o-po int ) , 

1  is  the  curlsquared  operator,  and  the  dimensionless 
variables  are  those  used  in  the  Vinghybrid"  code,  see 
Cornell  univ.  lps  rept.  no.  2SB  (fr iedman,  sudan,  denavit). 
note  that  it  is  possible  to  arrange  for  the  output  to  appear 
in  any  chosen  units  by  setting  bnot  appropr  iatel y .  the  relation 
is  athetaCouter  wall)  =  bnot  *  nr  /  2.  thus,  one  need  only 
figure  out  what  athetaCouter  wall)  is  in  one’s  units,  and 
then  compute  the  correct  bnot. 

in  these  equilibria  the  current  is  carried  entirely  by  hot  ions  with 
gyroradii  which  may  range  from  infinitesimal  (as  in  the  hill  vortex 
equilibria)  to  of  order  the  system  size  (as  in  ion  ring  equilibria), 
a  representat ion  of  the  distribution  function  is  obtained  by  assigning 
locations  and  velocities  to  particles  in  a  manner  which  yields  the 
correct  density  and  mean  azimuthal  velocity  in  each  grid  cell,  within 
small  errors  associated  with  finite  particle  size  and  the  limited  number 
of  particles  employed. 

code  usage  and  input  variables  are  described  on  comment  lines  in  the 
source  i t sel f . 


REFERENCES-the  f  irst  two  references  describe  the  program  and  the  normalization, 
respectively,  though  the  latter  is  not  difficult  to  figure  out  from 
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the  source,  and  the  code  is  largely  sel f -documented  with  comment  lines, 
the  latter  references  provide  useful  background  information  on  field 
reversed  equilibria. 


a.friedman,  proc.1930  Sherwood  meeting  on  theoretical  aspects  of 
controlled  thermonucl ear  research  (tucson,  arizona). 

a. friedman,  r.n. Sudan,  and  j.denavit,  "a  linearized  3d  hybrid 

code  for  stability  studies  of  field  reversed  ion  rings," 

cornel  1  university  laboratory  of  plasma  studies  report  no. 268  (1979, 

submitted  to  j . comp . phys . ) 

1. sparks,  j.m.finn,  and  r.n. Sudan,  bul 1 . am. phys . soc . 24,  955  (1979). 

b. marder  and  h.ueitzner,  plasma  phys. 12,  435  (1970). 

d . v. anderson,  j. Killeen,  and  rn. e . rens ink,  phys. fluids  15,  351  (1972). 
r . v. 1 ovel ace,  d . a. 1 arrabee,  and  h . h . f 1 e ischmann,  phys. fluids  22,  701 
(1979) . 

AVAILABILITY- 

users  can  obtain  a  copy  of  the  source  from  filem  directory  .takeme 
of  user  number  1234: 

filem  rds  1234  .takeme  alwith.  rr(esc)end  /tv 

the  file  thus  obtained  will  be  named  rrmmddyy,  where  mm  is  the  month, 
dd  is  the  day,  and  yy  is  the  year  that  the  source  was  created, 
the  source  can  be  compiled  using  chatr;  instructions  appear  on 
comment  cards  within  the  source.  to  list  the  source,  type: 

allout  hsp  <source>  ccsp.  seq.  box  <boxnumtoer>  rigidrotor  /tv 
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banner  <usc>  <source>  ccsp.  seq.  box  <boxnumber>  rigidrotor  /tv 

uhere  <usc>  i3  the  appropriate  user  service  center  designator, 

<source>  is  the  name  of  the  source  file  retrieved  from  filem,  and 
<boxnumber>  is  the  appropriate  box  number. 

DISTRIBUTION-  uni imited 


F.  Radial  Code  Notes  (R,R0 ,RZ,R0Z) 

C.  K.  Birdsall 

Method  C 

Radial  weighting  may  also  be  done  to  advantage  with  particles 
having  charge  density  which  is  not  uniform  within  the  particle.  Niels 
Otani  has  pointed  out  several  advantages  of  using  such  as: 

having  the  center  of  charge  and  mass  at  the  arithmetic  center;  producing 
the  familiar  symmetric  triangular  particle  shape  when  used  with  linear 
weighting  independent  of  r;  having  a  simple  weighting  algorithm  (less 
operations  than  in  Methods  A  and  B).  Let  us  see  these  in  detail. 


As  the  charge  per  unit  radius  is  p(r)27trdrdz  and  p(r)/v1/r,  this  quantity 
is  constant.  Because  charge/dr  is  flat  (as  it  would  be  in  an  x-y  grid), 
weighting  charge  to  the  grid  r . 1 s  can  be  done  as  in  x-y  grids.  That  is, 
the  charge  to  be  assigned  to  r.  of  particles  1,  by  linear  weighting,  is 


m 


where  f  is  the  fraction  of  the  nominal  particle  in  the  cell  about  r . , 


P“i 


P-i 


and  dV»rdrd4>dz.  So,  with  p  1/r 


f  =  (rjH’rp-i)/(rj+rrj-i)  "  a/Ar  *  {rj+rrP)/Ar 


Simi larly. 


q(r.  J  -  q  (b/Ar)  -  q  (r  -r  )/Ar  . 
j+1  p  P  P  J 

Note  that  there  is  no  division  by  r.  of  r.x1  or  r  as  was  true  in  Method 
A  and  B,  simplifying  the  weighting.  Note  also  that  the  particle  shape  as 
seen  by  the  grid,  that  is,  q(at  r.  as  r^  is  varied)  is  simply  the  same  as 
with  linear  weighting  to  rectangular  grids.  Thus,  what  is  known  about  S(x) 
particle  shapes  fits  $(r).  Also,  the  symmetric  digital  filtering  used  in 
rectangular  grids,  like  ( 1,2,1)  and  compensation,  like  (-1,6-1),  may  be  used 
di rectly . 

The  0  weighting  is  simply  linear  interpolation,  to  be  used  as 
already  given  in  Method  B. 

Hence,  in  2d,  full  r,0  grid.  Table  1  (of  QPR  IV,  1979)  is  sim¬ 
plified;  the  weights  are 

f  weighted  ar«»  |  _  farortar  -  ir)lq 
l  total  weighted  area  J  ArA<f> 

as  in  rectangular  weighting.  The  charge  weighting  holds  down  to  and  at  the 
origin,  with  no  special  consequences  .as  rp-*0  or  rp*®»  The  nominal  particle 
density  blows  up  as  1/r,  but  the  total  charge  of  the  rod  or  ring  does  not. 
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B.  I.  Cohen  notes  that  the  charge  density  weighting  of  Method  A, 
as  used  in  SUPERLAYER  etc.,  is  proper  for  use  with  Poisson's  equation,  where 
charge  density  is  wanted  at  the  grid  points,  and  not  q . .  (Finding  charge 
at  the  grid  points  is  more  for  use  with  Gauss'  Law  where  charge  inside  some 
volume  is  accumulated. )  Method  A  then  uses  linear  weighting  directly, 
weighting  particle  density  to  the  grid.  (The  exercise  given  in  the  write-up 
of  Method  A  on  calculating  the  charge  on  the  grid  serves  to  check  on  charge 
conservation. ) 

QPR  HI  1979,  p.  92  has  pp  -  q^/ (vol ume  of  charge)  then  gives 

p  *q  /27Tr  or6z.  This  volume  appears  to  be  valid  only  in  the  limit  6r  +  0. 

P  P  P 

The  exact  volume  is 


vo  I  ume 


/ 


outer 


i  nner 


if 


rdrd<(>dz 


2,tt<Sz 


r.)27rdz  . 

i 


If  the  particle  center  is  r  *  (r  +  r.)/2  and  particle  thickness  is  r  -r.*5r, 

p  o  i  o  i 

then  the  volumes  are  the  same,  with  no  need  to  consider  <5r-*0. 
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G.  POLARES:  A  Two-Dimens ional  Electrostatic  R-0  Code 
Niels  F.  Otani  (Prof.  C.  K.  Birdsall) 

While  it  is  often  convenient  and  acceptable  to  approximate  cylin¬ 
drical  systems  with  rectangular  models,  there  are  many  cylindrical  systems 
which  exhibit  effects  not  found  in  their  rectangular  counterparts.  For  ex¬ 
ample,  the  diocotron  instability  boundaries  for  a  non-neutral  hollow  cylinder 
of  charge  depend  explicitly  on  the  inner  and  outer  radius,^  Clearly  the  non¬ 
linear  evolution  of  such  a  system  would  also  be  of  some  interest,  and  is 
easily  studied  using  computer  simulation. 

With  this  in  mind,  we  are  in  the  process  of  developing  a  two- 
dimensional  code  which  is  based  completely  on  the  r-9  coordinate  system 
(polar-electrostatic  =  POLARES) ,  for  use  on  the  CRAY.  An  appropriate  weight¬ 
ing  scheme  has  been  used  to  accumulate  the  charge  density  and  fields  on  a 
radial  grid,  and,  partly  as  an  experiment,  a  grid  of  Fourier  amplitudes  has 
been  used  in  the  azimuthal  direction.  Both  guiding-center  and  ful 1 -dynamics 
time-centered  leapfrog  movers  have  been  incorporated  in  the  program  in  such 
a  way  as  to  make  them  interchangeable  for  either  species.  A  number  of 
graphics  packages  have  also  been  written.  All  are  easily  inserted  into  or 
removed  from  the  main  code  or  post-processor.  Some  allow  interactive 
graphics  viewing  via  DISSPLA.  Further  “modular"  diagnostics  packages  are 
contemplated. 

THE  CODE 

The  program  POLARES  may  be  thought  of  as  divided  into  five  sub¬ 
units:  a  main  source  (MAIN-S),  a  cl  iche  f  i  1e  ( INPCOMP) ,  a  binary  module 
library  (POLARLIB),  the  executable  file  (POLARES),  and  an  input  file 
(INPOLAR). 
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MAIN-S  is  a  generally  simple  Fortran  source  usually  consisting  of 
common  blocks  and  calls  to  various  packages  in  POLARLIB.  This  format  allows 
the  user  to  insert  various  movers,  charge  weighting  schemes,  background  charge 
distributions,  diagnostics,  graphics  packages,  etc.,  with  ease,  thus  contri¬ 
buting  to  the  flexibility  of  the  code.  One  might  suppose  that  this  format 
results  in  a  slow  code;  indeed,  we  find  that  the  running  speed  is  around 
500,000  particle-time-steps  per  minute  of  computer  time  using  a  full-dyn¬ 
amics  mover  with  24  radial  grid-points  and  24  azimuthal  Fourier  modes.  This 
is  not  particularly  fast,  but  we  have  not  made  full  use  of  the  vector! zat ion 
feature  of  the  CRAY  at  the  present  time.  Vector i zat ion  of  the  mover  and 
charge  accumulator  should  yield  significant  savings  in  time. 

Common  blocks  and  grid  dimensions  are  loaded  into  MAIN-S  and 
modules  of  the  binary  library  POLARLIB  during  a  precompiling  sweep  from 
information  supplied  by  the  cliche  file  INPC0MP.  Thus  the  number  of  parti¬ 
cles,  grid  points,  and  Fourier  modes  are  easy  to  change. 

POLARLIB  contains  all  the  main  building  blocks  of  the  simulation. 
These  include  graphics  packages  compatible  with  P0LARES,  diagnostics,  grid 
charge  accumulators,  movers,  initializing  routines,  a  particle  loader,  and 
a  field  solver.  Further  development,  refinement,  and  testing  of  most  of 
these  rout ines  is  envisioned;  additional  diagnositcs  and  graphics  packages 
are  also  expected  to  be  created. 

Some  of  the  more  important  routines  are  outlined  here. 

RH0EWT  and  RH0IWT  accumulate  charge  density  on  the  grid  from 
particle  electrons  and  ions  respectively.  One  may  envision  the  charge 
weighting  scheme  as  follows.  We  first  lay  on  the  polar  plane  a  series  of  con¬ 
centric  circles  each  A r  outside  the  one  before  it,  out  to  a  maximum  radius  a. 
These  represent  the  grid  points,  probably  more  properly  called  grid  circles. 
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Each  charged  particle  with  coordinates  (r.,0.)  is  separated  into  two  charged 
particles  of  charge  (rj+1 -r . /Ar)q  and  (r.-rVAr)q  located  at  (r^,0.)  and 
(rj+1,9.),  respectively.  Here  r .  and  r.+^  are  the  radii  of  the  two  nearest 
grid  circles,  and  q  is  the  charge  (per  unit  length),  (By  charge  particles 
in  this  code  we  of  course  refer  to  line  charges  extending  infintely  in  both 
directions  parallel  to  the  cylindrical  axis  of  symmetry.)  Each  fractional 
line  charge  is  now  considered  as  a  sum  of  sinusoidal  waves  existing  only  on 
its  own  grid  circle  according  to  the  formula 

a.  /  H  \ 

p(0)  »  27TrT  \  *  *  4-*  ^cos  cos  *  sin  m6j  s‘n  m9)l 

m- 1 

where  q1  is  the  fractional  charge  of  the  sub-particle  located  on  the  grid 
circle  of  radius  r1 .  The  weighting  factor  1 / r 1  arises  from  the  formula 
for  a  delta  function  in  polar  coordinates: 

CO 

62(x-x.)  *  — 5(r-r.)(1+2  y'  cos  m(9-0.))  . 

2lrfi  '  n*-1  ' 

The  9  shape  function  will  be  made  more  general,  essentially  just  changing 

2 

the  coefficients  in  the  sum,  e.g.,  to  drop  off  as  1/m  ,  implying  a  broader 
particle,  say,  A9  wide,  where  A9  is  something  like  Thus  Fourier 

amplitudes  of  the  charge  density  are  accumulated  for  each  of  the  mode  num¬ 
bers  up  to  M  for  each  of  the  a/Ar  grid  circles.  Notice  that  we  cannot  have 
a  grid  circle  of  radius  zero  in  the  scheme.  Therefore,  in  order  to  deal 
with  particles  inside  the  first  grid  circle,  we  have  devised  a  somewhat 
ad  hoc  scheme.  The  idea  is  to  replace  each  eligible  particle  with  an  equi¬ 
valent  charge  distribution  on  the  first  grid  circle.  By  "equivalent"  we 
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mean  that  the  new  charge  distribution  produces  the  same  potential  outside 
the  first  gri  circle  as  did  the  charge  it  replaces.  Thus  for  particles 


with  r.  < Ar  we  use  the  Fourier  amplitudes  implied  by  the  charge  distribu¬ 
tion 


M  r. 


-  (l  +  2  £  (”r”)  ^C0S  m6i  COS  ^  +  sinme^sin  me)^ 


This  scheme  has  the  advantages  that  (a),  the  long  range  force  is  correct 
and  (b),  a  particle  crossing  the  first  grid  circle  is  represented  by  con¬ 
tinuously  varying  charge  density  quantities.  Figure  (la)  illustrates  the 
charge  density  resulting  from  a  single  point  charge. 

Once  a  charge  density  has  been  established  on  the  grid,  the 
next  step  is  to  solve  for  the  potential  and  electric  fields.  This  is  accom¬ 
plished  by  the  subroutine  GRFIELDS.  The  task  of  GRFIELDS  is  fairly  straight¬ 
forward;  a  tridiagonal  solver  is  used  to  find  the  Fourier  components  of  the 
potential  from  the  finite  difference  version  of  the  Poisson  equation 
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In  the  last  formula,  the  c/s  notation  is  meant  to  indicate  that  the  cos  m0 


coefficients  of  (E^).  are  derived  from  sin  m0  coefficients  of  <t> .  ,  and 
9  j m  jm 

similarly  the  sin  m0  coefficients  are  derived  from  the  cos  m0  coef f i c i ents , 

with  a  change  of  sign.  The  multiplication  by  m  (for  V$  angular  component) 

may  be  replaced  by  something  like  mdif(mA0)  to  imply  a  more  local  gradient; 

also,  both  components  might  be  weighted  to  reduce  angle  errors  in  E;  these 

techniques  are  common  to  xy  codes. 

At  present  the  boundary  conditions  imposed  are  (a)  *0  at 

jm 

r»a  for  all  m  and  (2)  $ .  *0  at  r*0  for  all  m^O  and  <l>  *  <t>  for  m*0. 

jm  Qm  1m 

This  corresponds  to  a  system  surrounded  by  a  conducting  cylindrical  wall  of 
radius  a. 

The  electric  fields  and  potential  are  illustrated  in  Figure  1. 
Notice  the  ripple  in  and  E^  at  the  radius  of  the  point  charge.  This 

is  due  to  the  ripple  in  the  charge  density.  It  is  anticipated  that  smoothing 
in  the  Fourier  space  of  9 ,  or  using  a  line  charge  broadened  into  a  finite 
size  rod  will  eliminate  this  undesirable  effect. 

The  particles  are  advanced  to  their  new  positions  by  one  of  the 
movers  available  from  P0LARL1B.  Thus  far  electron  and  ion  guiding-center 
movers  MOVEGC  and  MOVIGC  and  ion  f ul l -dynami cs  mover  MOVIES  have  been  written. 
All  movers  find  the  electric  fields  at  each  particle  by  reversing  the  steps 


RHOIWT  and  RHOEWT  used  to  compute  the  charge  density.  That  is,  for  a  parti 
cle  located  at  (r.,0,),  the  electric  fields  are  first  found  at  (rj,9.)  and 
(rj+^,0.)  by  (for  5  function  weighting  in  9) 


and  then  are  interpolated  to  the  particle  by  linear  weighting 


E(r.  ,9.) 


The  guiding-center  mover  uses  this  value  of  E 
guiding-center  velocity  from  the  equation 


to  compute  a 


v  =  c-ExB  +  c^gxB 
'new  1  '  2  3 


where  c^  and  are  appropriate  constants. 

The  f ul 1 -dynami cs  mover  obtains  the  new  velocities  from  the 
old  velocities  using  a  method  analogous  to  the  one  employed  in  A.  B.  Lang- 
don’s  one-dimensional  electrostatic  code  ESI: 
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r  l  c  y  i  < 

(vr)  t+i-At 

■  <V2  *  f  <Er>,  T 

t+*At 

’  (V2  *  J(E»>t  T 

where  to  =  qB  /me. 
c  o 


Both  movers  then  move  particles  to  their  new  positions  accordii 


'(r  +(v  )  At)2  +  (vj2  At2 
t  r  new  9  new 


1  (v  )2  At2 

%  r,  +  (vp)  At  *  l  9  new 
t  r  new  2  r 


(v  )  At 

9  +  arctan  — - rr 

C  rt+tvr)newAt 


%  ^V9^newAt  L  +  ^vr^newAt 

rt  V 


where  v  = v  +At/2  for  the  ful 1 -dynami cs  mover, 
new  -  t 

The  approxi  mat  i  ons  given  are  valid  for  vAt/r  « 1  .  At  present  we 
use  the  approximate  formulas  for  particles  with  r^  >  Ar  and  rt+At>Ar*  If 
the  approximate  formula  yields  < Ar  or  If  r^  < Ar,  the  exact  formulas 


are  used. 


Finally  the  new  velocities  of  both  guiding-center  and  full-dyn¬ 


amics  movers  must  be  expressed  in  the  new  local  coordinate  system: 


v  +  v  cos  A9  +  vrt  sin  A9 
r  r  9 


va  -►  -v  sin  A9  +  v  cos  A8 
y  r  9 


where  A9  *  0  -0  .  Again,  away  from  the  origin  we  use 


cos  A9  %  1 


,  sin  A9  %  A9 
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The  offset  of  the  velocities  from  the  positions  in  the  full-dyn¬ 


amics  mover  is  necessary  for  the  mover  to  be  time-centered.  The  offset  is 
accomplished  at  the  beginning  of  the  run  by  the  routine  SETIV  for  ions  (and 
SETEV  for  electrons,  to  be  written),  by 


v 

ro 


_  a.  E  £t 

m  rO  2 


(Vi 


0o 


-  a  c 

m  60  2 


(v  )  .  =  (v  ) t  cos 

r  -At/2  r  1 


a)  At 

c 


(v9}1  Sin 


u)  At 

c 


(V-4t/2 


tu  At 


to  At 


=  (vr)1  Sin  +  (ve),  cos  -y 


Initial  particle  positions  and  velocities  are  read  into  the  code 


from  an  input  file  (INPOLAR)  by  the  P0LARL1B  module  LOADR. 
PRELIMINARY  RESULTS 


The  code  has  been  tried  on  a  few  systems,  the  time  evolutions 
of  which  are  known. 

Figure  (2)  shows  one  particle  moving  in  the  field  of  another 
fixed  particle.  There  is  also  a  uniform  external  magnetic  field  present. 
The  charge  does  not  feel  the  effects  of  its  own  wall  image.  We  expect  the 
ExB  motion  of  the  particle  to  be  along  an  equi potent i a  1  and  to  be  fastest 
where  |E|  is  largest.  The  computed  trajectories  are  consistent  with  these 
expectations.  Notice  the  effect  of  the  ripple  of  the  fixed  charge,  especi¬ 
ally  evident  when  the  gu id i ng-center  mover  is  used. 
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Figure  (3)  shows  some  snapshots  of  an  instability  in  a  ring  of 
.  cha/ge  in  a  uniform  magnetic  field,  simu  1  atad  h*s  i  ng  the  gu  i«d  i  ng-center  mover. 
Theory  for  the  diocotron  instability  predicts  all  modes  below  m  =  6  should  be 
unstable  for  our  parameters.^  This  seems  to  be  consistent  with  our  simula¬ 
tion. 

Figure  (4)  illustrates  the  charged-sheet  instability  which  has 

2 

been  observed  in  electron  beams  and  is  believed  to  occur  in  auroral  displays. 
The  mover  used  here  is  the  f u 1 1 -dynami cs  mover,  but  the  guiding-center  mover 
produced  qualitatively  very  similar  results.  This  allows  us  to  be  somewhat 
confident  that  both  movers  are  working  correctly. 

Figure  (5)  illustrates  our  attempt  to  simulate  the  Rayleigh- 
Taylor  instability  of  a  neutral  plasma  moving  in  a  gravi tat ional  field  di¬ 
rected  radially  outward,  g  simulates  either  VB  or  R  drifts,  which  are  signed, 
hence  lead  to  flute  instabilities.  The  ions  are  pushed  with  a  ful 1 -dynami cs 
mover  while  a  guiding-center  mover  is  used  for  the  electrons.  We  observe  some 
peculiar  behavior  in  the  vicinity  of  the  origin.  This  is  quite  likely  an 
artifact  of  the  code  and  will  be  investigated  in  the  near  future.  As  is 
obvious  in  the  figure,  the  growth  of  the  instability  is  quite  rapid.  When 
gravity  of  the  same  magnitude  is  directed  inward,  the  plasma  is  found  to  be 
stable.  We  also  find  a  mild  growth  rate  even  in  the  absence  of  gravity.  This 
is  possibly  due  to  expansion  of  the  warm  plasma  (in  uniform  B)  or  the  abnor¬ 
mally  large  field  fluctuations  likely  to  be  associated  with  so  few  particles 
(1024  ions,  1024  electrons).  The  large  time  step  used  in  this  preliminary 
run  (u>  At*l)  is  also  a  possible  cause. 

Finally,  plasma  oscillations  have  been  observed  in  a 
cylindrical  plasma.  In  this  run,  the  electrons  were  pushed  using 
the  f u 1 1 -dynam i cs  mover  with  no  external  magnetic  field  present,  and 
the  ions  were  fixed.  The  frequency  of  oscillation  was  not  correct  however; 
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we  found  oscillations  occurring  at  about  ui  ,  instead  of  the  theoretical 

Pe 

value  of  w  //2,  the  dipole  resonance  frequency.  (Actually,  in  the  presence 
pe 

of  the  wall,  and  with  our  parameters,  the  theoretical  resonant  frequency 
is  lowered  slightly  to  No  explanation  for  the  discrepancy  is 

offered  at  this  time. 


FUTURE  DEVELOPMENTS 


The  code  POLARES  has  thus  far  produced  encouraging  results 
though  most  are  as  yet  qualitative  in  nature.  It  is  hoped  that  further 
work  with  the  code  will  produce  quantitative  results  and  clear  up  problems 
already  encountered. 

Definitely  on  the  agenda  are  vector i zat i on ,  which  should  result 
in  considerable  savings  in  computer  time,  smoothing  in  9-Fourier  space, 
and  improvement  in  the  diagnostics.  Results  of  these  developments  will 
be  reported  in  a  future  QPR. 
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FIGURE  CAPTIOUS 


FIG.  1. 


FIG.  2. 


(a) 


(b) 

FIG.  3. 


FIG.  4. 


FIG.  5. 


Charge  density,  potential,  and  electric  fields  (in  arbitrary  units) 
generated  by  a  line  charge  located  at  r*l6Ar,  0*tt/2.  The  radial 
grid  appearing  in  the  plots  is  real,  but  the  9-grid  is  generated 
by  the  plotting  routine,  a  a  2^Ar,  M  =  24. 

Orbit  of  a  positively-charged  particle  in  the  field  of  a  fixed  posi¬ 
tively-charged  particle  lor^Led  at  r»l6Ar,  9  *  tt/2.  a  =  24Ar,  M  =  24. 

2  2 

Full-dynamics  mover:  In  cgs,  4irq  At  /m=2,  qB  At/mc*  0.5,  1500  time- 

o 

steps.  Initial  conditions  for  the  mobile  particle:  r»8Ar,  9»tt/2, 

v  At  =  0,  v.At  =  0.667. 
r  d 

Gu  i  di  ng-center  mover  :  in  cgs,  47rqcAt/BQ  *  2  ,  3000  timesteps.  Initial 
conditions:  r  =  6.667Ar,  0=tt/2. 

Instability  in  a  charged  ring  using  the  guiding-center  mover,  1024 
particles  uniformly  distributed  between  r  *  9Ar  and  r*1lAr,  4nqcAt/BQ 
=  0.02,  a  *  l6Ar,  M  *  16. 

Charged-sheet  instability  using  the  cu 1 1 -dynami cs  mover  with  1024 

particles,  uj  At  38  0.9,  At  =  1  ,  a  =  24Ar ,  M  *  24. 

P  c 

Simulation  of  a  un i forml y-magnet i zed  plasma  in  a  gravitational  field 

directed  radially  outward.  A  full  dynamics  ion  mover  and  an  ExB 

guiding-center  electron  mover  are  used.  Initially  the  plasma  is 

uniformly  distributed  at  radius  8Ar.  Parameters:  w  . At  ■  1 , 

2  2 

up ;  -0.22,  wpe/|ci>ce  |  ■  “p|/wcf .  g  ■  0.0625(r/Ar)  grid  points/(At)^, 

1024  electrons,  1024  ions,  a*24Ar,  M*24. 
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■H.  •  •  Dispersion  Re  lat  i  ott*  9o*l  ver  For  •a**L*rnreari  zeth-Two-Fl  u  i  &  Urri  fcrrm  Pl-asma 
Model 

Niels  F.  Otani  (Prof.  C.  K.  Birdsall) 

A  simple  solver  is  being  developed  for  the  dispersion  relation 
associated  with  a  linearized,  magnetized,  electromagnetic,  two-fluid,  uni¬ 
form  plasma  model.  The  model  allows  an  isotropic  temperature  for  each  species. 
The  unperturbed  quantities  nQ,  T  ,  and  Bq  are  uniform  in  space  and  constant 
in  time,  and  neither  species  is  allowed  to  drift.  The  linearized  equations 
required  to  describe  this  system  are  then  as  follows: 


+  n  7  •  5v 
o  -s 
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where  s  refers  to  the  species  involved  (here  electrons  and  ions).  This  code 
is  intended  as  an  educational  aid  only;  the  zoo  of  sma 1 1 -ampl i tude  waves  and 
oscillations  resulting  from  this  set  of  equations  has  been  wel 1 -understood 
for  some  time. 

Allowing  the  perturbed  quantities  to  have  the  dependence 
exp ( i (k« r  •  wt) ) ,  the  dispersion  relation  for  this  set  of  equations  is  easily 
found  to  be* 
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The  dispersion  relation  may  be  expressed  as  a  30th-degree  poly¬ 
nomial  in  w,  but  generally  most  of  the  roots  are  zero.  The  code  therefore 

samples  the  value  of  the  polynomial  at  ur  ^  amd  2u>  ^  where  u>  is  much 

test  test  test 
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smaller  than  any  frequency  intrinsic  to  the  free  parameters  of  the  system. 


Then  log?  f(2w  ) /f7 (u>  )  factors  of  u)  are  divided  out  of  the  polyno- 

4  L v j L  LCD L 

mial  f,  and  the  new  function  is  used  in  Steven  Au-Yeung  and  Alex  Friedman's 

2 

root  solver  code  SOLVER  to  obtain  the  various  branches  of  w(|k|)  (9  fixed). 
This  procedure  is  necessary  as  SOLVER  often  has  trouble  with  functions  hav¬ 
ing  several  equal  roots.  An  example  of  the  output  is  shown  in  Fig.  1. 

Future  development  of  the  code  will  allow  plotting  of  oj(9) 

( | k |  fixed)  and  also  permit  the  user  to  obtain  relative  amplitudes  and 
phases  of  5E,  5B,  5 n,  6v,  etc.,  which  provide  information  about  the  nature 
of  the  mode.  Either  the  subroutines  used  with  SOLVER  or  a  complete  code 
will  be  made  available  to  users. 


REFERENCES 


1.  Allis,  W.  P.,  Buchsbaum,  S.  J.,  and  Bers,  A.,  Waves  in  Anisotropic 
PI asmas .  Ch .  1.  M.l.T.  Press,  1 96 3 • 

2.  Au-Yeung,  H.  Stephen,  and  Friedman,  Alex.  "SOLVER:  An  Analytic 
Function  Root  Solving  and  Plotting  Package."  Electronics  Research 
Laboratory:  University  of  California,  Berkeley,  August,  1979*  Memor¬ 
andum  No.  UCB/ERL  M79/55. 


I.  RJET  Development 


No  special  report. 
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Software  Developments 
Niels  F.  Otani 


A  utility  routine.  DD38MAKE.,  has  been  created  for  use  cn  the  A-machine 
for  the  purpose  of  facilitating  the  conversion  of  frS3  files  into  ddoO 
files.  The  code  is  especially  useful  when  several  f r 30  files  have  been 
generated  on  the  A-machine  or  sent  from  the  Cray  to  be  viewed  as  dd89  files 
using  72i<PL0T.  The  routine  may  be  obtained  from  F  IL£;I  by  typing: 
filem  rds  1213  . public  dd83mahe 

DDSQihAKE  is  a  L!2-creat©d  1  ibrar'.;  or:!  contains  the  executable  file  c.s 

t  re  f  irst  entry.  Also  included  in  the  1  ibra;**.:  are  the  scurce  and  the  utility 

routine  DDcX.  The  latter  is  a  ddS0  file  editor  described  in  a  previous 
1 

UPR  which  is  often  useful  in  conjunction  with  this  program. 

The  execute  line  for  DD30MAKE  is: 
dd30mal<e  box  bnn  any  id  l  end  ]  /  t  v 

If  the  box  and  id  are  omitted..  DP30MAKE  will  prompt  for  this  information. 

If  11  end"  is  omitted  in  the  execute  line.  c!c!?9ma.'<©  will  attempt  to 
convert  any  files  beginning  with  "  f  13^  '  into  :JdS9  files.  The  new  ddE3 
files  will  retain  the  two  identifying  characters  c f  their  old  frSG 
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Station  III 

PLASMA  SIMULATION  TEXT 


A  contract  has  been  signed  with  McGraw-Hill  for  the  text  to 
be  delivered  to  them  Dec.  1,  1980,  duly  shortened  and  updated.  No 
further  progress  will  be  reported  here. 

Station  IV 

SUMMARY  of  REPORTS,  TALKS,  PUBLICATIONS 

Abstracts  of  papers  to  be  presented  at  the  Sherwood  Fusion 
Theory  Meeting  April  2 3 “ 25  at  the  University  of  Arizona,  Tucson  follow. 

A  listing  of  papers  dealing  with  instabilities  due  to  cur¬ 
rents  in  diode  regions  is  also  given. 
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FRM  and  ION  RING  EQUILIBRIA  for  the  RINGHYBRIO  CODE* 

Alex  Friedman 

Electronics  Research  Laboratory 
University  of  California 
Berkeley,  CA  9^720 

We  have  developed  a  new  method  of  generating  axi symmet r i c  field 
reversed  equilibria  for  the  RINGHYBRID  code,  a  linearized  3d  hybrid  code 
used  for  studies  of  the  1 ow-f requency  stability  of  such  configurations1*2. 
Equilibria  of  known  functional  form  are  calculated  numerically,  and 
appropriate  particle  initial  condition  and  magnetic  field  data  are  used 
to  initialize  the  hybrid  code. 

Exponential  rigid  rotor  Vlasov  equilibria,  with  distribution  functions 
of  the  form  f « exp [- (H-QPg)/T],  are  computed  on  an  r-z  grid  using  an 
iteration  scheme  similar  to  that  of  Sparks  and  Finn3.  In  these  equilibria 
the  current  is  carried  entirely  by  hot  ions  with  gyroradii  which  may  range 
from  infinitesimal  (as  in  the  Hill  vortex  equilibria)  to  the  order  of  the 
system  size  (as  in  ion  ring  equilibria).  A  representation  of  the  distribu¬ 
tion  function  is  obtained  by  assigning  locations  and  velocities  to  particles 
in  a  manner  which  yields  the  correct  density  and  mean  azimuthal  velocity  in 
each  grid  cell,  within  small  errors  associated  with  finite  particle  size  and 
the  limited  number  of  particles  employed  (several  thousand).  RINGHYBRID 
advances  these  particles  along  their  equilibrium  orbits  in  the  time- 
independent,  self-consistent  zero  order  field.  (For  stability  studies  the 
code  simultaneously  advances  linearized  displacements  from  these  orbits,  as 
well  as  linearized  quasineutral  fluid  equations  for  cold  ion  and  electron  com¬ 
ponents,  and  self-consistent  first  order  E  and  B  fields,  which  are  fully  3d*) 

A  variety  of  field  reversed  mirror  and  ion  ring  equilibria  have  been 
obtained;  when  the  equilibrium  field  and  particle  initial  condition  data  are 
used  in  RINGHYBRID,  particle  moments  remain  nearly  constant  in  time.  Orbits 
are,  in  general,  quite  complicated,  and  preliminary  indications  are  that 
many  are  truly  stochastic2.  The  complexity  of  modera te-gy rorad i us  FRM 
orbits  arises  from  the  fact  that  particles  moving  along  closed  field  lines 
may  or  may  not  reflect  in  the  high  field  regions,  and  furthermore  may  pass 
near  the  field  null  during  their  gyratory  excursions.  Also,  the  rate  of  a 
particle's  azimuthal  drifting  motion  depends  strongly  upon  where  it  lies  on 
its  orbit  in  the  r-z  plane. 

Previously,  equilibria  were  generated  by  injecting  bursts  of  particles 
into  an  external  magnetic  field  and  solving  Ampere's  law  (with  a  resistive 
source  term  'tf03Ag/3t  included  to  damp  collective  oscillations)  repeatedly 
until  Ag  became  nearly  stationary.  While  this  produced  equilibria  which  were 
useful  for  ion  ring  stability  studies2,  the  new  method  affords  greater  control 
over  the  resulting  equilibria,  and  so  parameters  may  more  readily  be  varied. 

The  package  which  computes  equilibrium  fields  and  particle  initial  con¬ 
ditions  is  called  RIGIDR0T0R,  and  is  available  for  genera)  use  through  L1BRIS 
on  the  NMFECC  CDC-7600. 

A 

Research  supported  by  the  U.S.  Department  of  Energy  under  Contract  No.  DE-AS03- 
76SF0003^,  Project  Agreement  No.  DE-AT03-76ET53064 . 

lA.  Friedman,  R.  N.  Sudan,  and  J.  Denavit,  Proc.  Eighth  Conf.  on  Numerical 
Simulation  of  Plasmas,  Paper  No.  PC-13,  Lawrence  Livermore  Lab.  Conf.  Proc. 

No.  CONF-78O61I4  (1978);  Cornell  Univ,  Lab.  of  Plasma  Studies  Rept.  No.  268 
(1979,  to  appear  in  J.  Comp.  Phys.). 

2A.  Friedman,  J.  Denavit,  and  R.  N.  Sudan,  Bull.  Am.  Phys.  Soc.  956  (1979). 

3L.  Sparks,  J.  M.  Finn,  and  R.  N.  Sudan,  Bull.  Am.  Phys.  Soc.  2A,  955  0979). 
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SATURATION  OF  THE  LOWER-HYBRID  DRIFT  INSTABILITY* 


Yu-Jiuan  Chen  &  C.  K.  Birdsall 
Electronics  Research  Laboratory 
University  of  California 
Berkeley,  CA  94720 

Bruce  I .  Cohen 

Lawrence  Livermore  Laboratory 
University  of  California 
Livermore,  CA  94550 

Saturation  mechanisms  of  the  lower-hybrid  drift  instability  in  the  low 
drift  velocity  regime  (vd«  v  therma  1  i  n^  are  stuc^ec*  both  *n  a  Id  particle- 
hybrid  simulation  and  a  non l i nearaper?urbat i on  theory;  v,  is  the  relative 
drift  velocity  between  the  electrons  and  ions.  We  find  that  ion  traping  sta¬ 
bilizes  the  instability  when  v^  is  kept  constant,  unless  current  relaxation 
(vd+0)  occurs  via  quasi  linear  diffusion.  We  also  find  a  new  saturation 
mechanism:  a  nonlinear  ion  orbit  perturbation  induces  a  frequency  shift 
which  is  particularly  effective  in  saturating  long  wavelength  modes. 

We  use  a  one-dimensional  slab  configuration.  In  zero  order,  the  ion 
pressure  gradient  cancels  the  equilibrium  electric  field  force.  in  first 
order,  the  ions  are  treated  as  unmagnetized,  and  the  electrons  are  assumed 
to  respond  to  the  wave  linearly,  because  wcj  <*:  |  w  |  ^pe  *  wce . 

In  simulation  of  small  amplitude  modes,  there  was  good  agreement  of 
the  linear  growth  rate,  real  frequency,  and  influence  of  finite  beta  effects 
associated  with  the  nonresonant  gradient  B  electron  orbit  modifications 
with  linear  theory.  For  large  amplitude  modes,  at  zero  plasma  beta  and  zero 
electron  temperature,  the  simulations  show  that,  when  v^  is  kept  constant, 
the  lower-hybrid  drift  instability  is  stabilized  by  ion  trapping.  When  v<j 
is  allowed  to  vary  in  time  (for  example,  as  a  self-consistent  consequence 
of  momentum  conservation),  stabilization  occurs  as  result  of  current  relax¬ 
ation  with  v^  +  0,  at  a  much  lower  level,  a  little  below  Davidson’s  prediction 
for  v^/v^j  <  1  ,  as  given  by1 

i  /iV 

8(1  +  0}^  /(i>^  )  mi  \Vti  / 
pe  ce  \  / 

Analytic  theory  shows  that  a  finite  perturbation  of  the  ion  orbits  leads 
to  a  nonlinear  frequency  shift  that  can  stabilize  the  lower-hybrid  drift  in¬ 
stability.  However,  the  resulting  saturation  level  is  small  only  for  modes 
with  wavelengths  much  longer  than  that  of  the  most  unstable  mode.  This  re¬ 
sult  is  obtained  from  a  self-consistent  solution  of  the  VI asov-Poi sson 
equations  using  perturbation  theory  in  which  the  nonlinear  dielectric  func¬ 
tion  and  the  nonlinear  temporal  evolution  of  a  single  unstable  mode  in  the 
low  drift  velocity  regime  are  calculated  analytically. 

This  research  was  supported  in  part  by  the  Office  of  Naval  Research  under 
Contract  N000 1 4- 77~c-0578 ,  and  in  part  by  the  Department  of  Energy  under 
Contract  No.  W-7405-ENG-48 . 

!R.  C.  Davidson,  Phys.  Fluids  21,  1375  (1978). 
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INSTABILITIES  DUE  TO  CURRENTS  IN  DIODE  REGIONS 
(C .  K.  Bi rdsa  1 1 ) 

Interest  in  these  instabilities  has  varied  in  time;  some  claim 
that  this  interest  revives  on  about  a  seven  year  cycle,  beginning  with  the 
work  on  limiting  current  (or  perveance)  by  B.  Salzberg  and  A.  V.  Haeff  in 
the  late  1930's  and  that  of  J.  R.  Pierce  in  the  early  19^0's.  We  may  be 
near  another  peak  now.  While  later  cycles  have  tended  to  progress  beyond 
earlier  work,  there  is  evidence  that  later  workers  have  tended  to  miss 
earlier  work. 

Starting  in  1959,  our  plasma  theory  and  simulation  group  contri¬ 
buted  some  insights  into  these  problems.  One  insight  from  simulation  was 
that  the  diode  current  instabilities  could  grow  into  large  amplitude  oscil¬ 
latory  steady  states  Another  was  that  the  time-average  of  the 

oscillatory  current  was  larger  than  that  obtained  from  time  independent 
ana  lysis. 

A  list  follows  of  representative  (but  not  exhaustive)  papers 
of  our  group  and  related  work  done  elsewhere  in  the  1960's.  The  papers 
themselves  provide  more  references.  We  commend  them  to  those  beginning  to 
work  on  current  instabilities:  in  diode  regions  as  such,  or  in  sheaths, 
divertors,  guns,  accelerators,  propulsion  devices,  direct  converters,  etc., 
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InAtabilitieA  due  to  CuAAenl6  in  Diode.  Region  •  Selected  Ze{eA.ence6  in  1960*  6 


i960 

Lomax,  R.  J.,  Transient  space-charge  flow,  J.  Elec,  and  Cont.  9_,  pp.  127- 
140,  August. 


1961 

Birdsall,  C.  K.  and  Bridges,  W.  B. ,  Space-charge  instabilities  in  electron 
diodes  and  plasma  converters,  J.  Appl .  Phys.  34.,  pp.  2611  - 26 1 8 ,  December. 

Lomax,  R.  J.,  Unstable  electron  flow  in  a  diode,  Proc.  I.E.E.  Part  C,  108, 
pp.  119*121,  March. 

1962 

Bridges,  W.  B.  and  Birdsall,  C.  K. ,  An  electron  stream  instability, 

185  pages,  Tech.  Rept.  No.  60-443,  E.R.L.,  Univ.  of  Calif.,  Berkeley,  CA, 

March.  (Bridges'  Thesis) 

1963 

Bridges,  W.  B.  and  Birdsall,  C.  K.  ,  Space-charge  instabilities  in  electron 
diodes,  II,  J.  Appl.  Phys.  34,  pp.  2946-2955,  October. 

Buneman,  0.  and  Kooyers,  G.  P . ,  Computer  simulation  of  the  electron  mixing 
mechanism  in  ion  propulsion,  A.I.A.A.J.,  J_,  pp.  2525-2528. 

1964 

Burger,  P. ,  Nonexistence  of  dc  states  in  low-pressure  thermionic  converters, 

J.  Appl.  Phys.  35.,  pp.  3048-3049,  October. 

Cutler,  W.  H. ,  High  frequency  oscillations  in  a  thermal  plasma,  J.  Appl. 

Phys.  35_,  PP*  464-465,  February. 

1965 

Bridges,  W.  B. ,  Frey,  J.  I.  and  Birdsall,  C.  K. ,  Limiting  stable  currents  in 
bounded  electron  and  ion  streams,  IEEE  Trans.  Electron  Devices  J_2,  pp.  264-272, 
May. 

Burger,  P. ,  Theory  of  large  amplitude  oscillations  in  the  one-dimensional 
low  pressure  cesium  thermionic  converter,  J.  Appl.  Phys.  36^,  PP-  1938-1943, 
June. 

Burger,  P. ,  Dunn,  D.  A.  and  Halsted,  A.  $.,  Computer  experiments  on  the  ran¬ 
domization  of  electrons  in  a  col  1  i  s  ion  less  plasma,  Phys.  Fluids  8^,  pp.  2263- 
2272,  December. 

Frey,  J.  I.  and  Birdsall,  C.  K. ,  Electron  stream  diode  instabilities  with 
elastic  collisions,  J.  Appl.  Phys.  36_,  pp.  2962-2964,  September.  (See  cor¬ 
rection  by  Faulkner  and  Ware,  1969-) 

Wadhwa,  R.  P.,  Buneman,  0.  and  Brauch,  0.  F. ,  Two-dimensional  computer  exper¬ 
iments  on  ion-beam  neutra 1 i zation,  A.I.A.A.J.  1076-1081,  June. 
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1966 


Birdsall,  C.  K.  and  Bridges,  W.  B.,  Electron  Dynamics  of  Diode  Regions, 
Academic  Press,  N.  Y.  (See  Chap.  3*  Stability  of  flow;  nonlinear  solu¬ 
tion  of  mu  1 1 i part i cl e  model.). 

Brauch,  D.  F.,  Buneman,  0.  and  Wadhwa ,  R.  P.,  Computer  studies  of  plasma  boun¬ 
daries  and  lens  effects  created  by  immersed  and  withdrawn  neutralizers, 

A. 1 . A. A. J .  A,  651-653,  April . 

Cutler,  W.  H.  and  Burger,  P. ,  Oscillations  in  the  thermal  cesium  plasma 
diode,  J.  Appl .  Phys.  37,  PP-  2867-2873,  June. 

Frey,  J.  I.  and  Birdsall,  C.  K. ,  Instabilities  in  a  neutralized  electron 
stream  in  a  finite  length  drift  tube,  J.  Appl.  Phys.  36^  PP-  2962-2964, 

April  (See  correction  by  Faulkner  and  Ware,  1969.). 

Twombly,  J.  C. ,  Dynamic  behavior  of  a  long  thin  electron  beam,  IEEE  Trans,  on 
Elec.  Dev.,  ED-13,  PP-  934-942,  December. 

1967 

Burger,  P. ,  Elastic  collisions  in  simulating  one-dimensional  plasma  diodes 
on  the  computer,  Phys.  Fluids  10 ,  pp.  658-666,  March. 

Hirsch,  R.  L. ,  I nert ia 1 -electros tat i c  confinement  of  ionized  fusion  gases, 

J.  Appl.  Phys.  38.’  PP*  4522-4534,  October. 

1968 


Hockney,  R.  W. ,  Formation  and  stability  of  virtual  electrodes  in  a  cylinder, 
J.  Appl.  Phys  39_,  pp.  4166-4170,  August. 

1969 

Faulkner,  J.  E.  and  Ware,  A.  A.,  The  effect  of  finite  ion  mass  on  the  stab¬ 
ility  of  a  space-charge-neutral i zed  electron  beam,  J.  Appl.  Phys.  40,  pp. 
366-368,  January. 

1970 

Barnes,  C.  W. ,  The  computer  simulation  of  a  spherically  symmetric  plasma, 
SUIPR  Report  No.  344,  Institute  for  Plasma  Research,  Stanford  University, 
Stanford,  CA,  March. 

1971 


Dunn,  Donald  A.,  Models  of  Particles  and  Moving  Media,  Academic  Press,  New 
York  (See  Chap.  6,  Motion  of  many  interacting  particles  -  the  Lagrangian 
mode  1 . ) . 
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